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HIGHER-ORDER NUMERICAL SOLUTIONS USING CUBIC SPLINES 


S. G, Rubin and P. K, Khosla 

Polytechnic Institute of New York 
Farmingdale, New York 

ABSTRACT 

A cubic spline collocation procedure has recently been 
developed for the numerical solution of partial differential 
equations. In the present paper, this spline procedure is 
reformulated so that the accuracy of the second-derivative 
approximation is improved and parallels that previously ob- 
tained for lower derivative terms. The final result is a 
numerical procedure having overall third-order accuracy for 
a non-uniform mesh and overall fourth-order accuracy for a 
uniform mesh. Solutions using both spline procedures, as 
well as three-point finite difference methods, will be pre- 
sented for several model problems. 


I . INTRODUCTION 

1 2 

In a recent study Rubin and Graves ' have presented a 
3 4 

cubic spline ' collocation procedure for the numerical solu- 
tion of partial differential equations. This technique ex- 
hibits the following desirable features: (1) The governing 

matrix system is always tridiagonal so that well-developed 
and highly efficient inversion algorithms are applicable; 

(2) cubic spline interpolation leads to second order accuracy 


for second derivatives, e.g., diffusion terms in the Navier- 
Stokes equations. This order of accuracy is maintained even 
with rather large non-uniformities in mesh width; (3) first 
derivatives or convection effects are fourth-order accurate 
for a uniform mesh and third-order with mesh non-uniformity; 

(4) derivative boundary conditions can in many cases be applied 
more accurately and with less difficulty than with conven- 
tional finite-difference schemes; (5) a simple two-point 
relationship exists between the spline approximation for the 
first and second derivatives; and (6) unlike finite-element 
or other Galerkin (integral) methods, which are generally not 
tridiagonal, the evaluation of large numbers of quadratures is 
unnecessary. 

1 2 

Solutions have been obtained for a number of problems ' 
with explicit, implicit and spline alternating direction 
implicit (SADI) temporal or spatial marching procedures. 

Moreover, for the viscous and potential flow problems consid- 
ered, it was found that with the spline procedure there was 
no particular advantage gained with the equations in divergence 
form. In some recent studies it has been found that the 
divergence form may be desirable with flux boundary conditions. 
These results are described later in this paper. 

Agreement of the spline solutions with exact analytic 

results and very accurate finite-difference solutions obtained 

12 12 

with a very fine mesh has been quite good ' . All comparisons ' 
with conventional three-point finite difference formulations 


2 



demonstrate the improved spline accuracy associated with 
(i) the higher-order convection approximation, (ii) the treat- 
ment of derivative boundary conditions, or (iii) the higher- 
order accuracy of spline second derivatives (diffusion) when 
specifying a non-uniform mesh. Solutions for the Burgers 
equation, the two-dimensional diffusion equation and the in- 
compressible viscous flow in a driven cavity are found in 
Refs. 1 and 2. 

In the present paper, the cubic spline procedure is 
reformulated so that the accuracy of the second-derivative 
approximation is improved and parallels that obtained for the 
lower derivative terms. The final result is a combined spline- 
finite difference numerical procedure having overall third- 
order spatial accuracy for non-uniform meshes and overall 
fourth-order spatial accuracy with a uniform mesh. In order 
to differentiate the two spline procedures, we shall designate 
the original spline formulation ' as spline 2 and the improved 
formulation presented here as spline 4. 

As shown in sections II and III, the cubic spline colloca- 
tion procedure involves a third-order interpolation polynomial 
with the function and the second (or first) derivative of the 
function as unknowns at each mesh point. continuity of the 
first (or second) derivative leads to the tridiagonal system 
of equations to be considered. In section IV, it is shown how 
the familiar central difference second-order accurate finite- 
difference theory results from a quadratic spline interpolation 
procedure. Using the earlier spline designation, the finite- 
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difference theory is classified as spline 1. 

Recently, several higher-order finite-difference schemes 
with similar properties have been proposed; i.e., the functions 
and derivatives are considered unknown at each mesh point, or 
the functions are collocated at three points instead of one. 

The methods which have been termed Hermitian finite-difference^'^ 
Fade approximation"^ or compact differencing^, and Mehrstellung^ 
have been developed for a uniform mesh and have somewhat lower 
truncation errors than the five-point pentadiagonal fourth-order 
finite difference procedure. As with the spline formulation, 
they remain of tridiagonal form. 

The authors have examined these procedures, as well as a 
fourth-order spline-on-spline method, and found them to be, in 
fact, identical ; i,e., any one can be derived from any of the 
others. As with spline 4, these finite-difference or spline- 
on methods are fourth-order with a uniform mesh and third- 
order with a non-uniform mesh. The main differences are 
handling of the boundary conditions, the relationship between 
the approximations for the convection and diffusion terms, and 
the truncation errors. The truncation errors for first deriv- 
atives are identical. The truncation error for the second- 
derivative to be discussed later for spline 4 is 50% smaller 
than that found with the higher-order finite-difference or 
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spline-on- spline collocation formulae. 

In order to evaluate the spline procedures, the trunca- 
tion errors, stability limitations and effects of boundary 
conditions will be discussed. Spline 2 is reviewed in section 
II, spline 4 is introduced and discussed in section III, and 
spline 1 is presented in section IV. The stability conditions 
for all methods are outlined in section V. Solutions using 
both spline procedures, as well as a three-point finite-difference 
method, are presented for several model problems in section VI. 

Both uniform and non-uniform meshes are considered. In each 
case the analytic solution or a very accurate numerical solu- 
tion is available for comparison purposes. The problems to be 
considered include (1) a boundary layer-like solution of Laplace's 
equation, where a spline relaxation method is applied, (2) poten- 
tial flow over a circular cylinder with a spline successive 
approximation procedure, (3) the weak shock solution for the 
nonlinear Burgers equation by a two-step explicit or an implicit 
spline integration, (4) divergence and non-divergence solutions 
for the linear Burgers equation with flux and other derivative 
boundary conditions, (5) the impulsive motion of right angle 
corner (Rayleigh problem) with SADI, (6) the solution of the 
two-point boundary value problem describing similar boundary 
layer behavior, and (7) non-similar constant pressure boundary 
layer solutions for large Reynolds number using physical variables. 
The results will be summarized in section VII. 
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II. SPIiINE 2 - REVIEW OF CUBIC SPLINE THEORY 


Consider a mesh with nodal points such that 


a=x^<Xi <X2 




and with 


h . =x . -X . , > 0 

1 X x-1 


Consider a function u (x) such that at the mesh points x^^, 
u(Xj^)=u^. The cubic spline is a function S ^(u^, x) =S ^ (x) which 
is continuous together with its first and second derivatives on 
the interval [a,b], corresponds to a cubic polynomial in each 
sub-interval satisfies S ^(Uj^;x^) =u^ . In the 

usual spline terminology, spline 2 is defined as a cubic spline 
of deficiency one, since all but one of the three polynomial 
derivatives are continuous. 

If u(x) and its derivatives are continuous, it has been 
shown that the spline function S^(x) approximates u(x) at all, 
points in [a,b] to fourth order in max h^ . First and second 
derivatives of S^(x) approximate u'(x) and u"(x) to third and 

3 

second order, respectxvely . See Ahlberg, Nilson and Walsh for 
detailed proofs of convergence. 

If S^(x) is cubic on [x^_j^,x^], then in general. 




x-x , 


where M . = S (x . ) • 
X A x 


Integrating twice leads to the interpolation formula on 
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(x.-x)^ (x-x. , 

S,(x)=M. , + M.- 7 -t ‘— 

A x-1 6h. 1 6h. 


M. ^h? (x.-x) M.h? (x-x. ,) 

(“i-x- -t— -p) — h— 

1 1 


(la) 


The constants of integration have been evaluated from 

and S.(x. ,)=u. S . (x) on [x.,x.,,] is obtained with i+1 

A 1-1 1-1 A 1 1+1 

replacing i in (la) , 

The unknown derivatives are related by enforcing the 
continuity condition on S^(x). With S^(x7)=m7 on [i-l,i] and 
S^(x^)=mi ^^i'^i+1^' require m” = mt = m^ . We find for 

i“l # 


h. h.+h. T h. , 

■— M. + M.+ — M. ,T- , 

6 1-1 3 1 6 1+1 h 




i+1 


(u . -u . T ) 
1 i-l'^ 

h. 

1 


(lb) 


Additional spline relationships that are easily derived are 
listed below; 


^ m^. ^t+ 2(^^ + tn . + — m. 


h. “i-r-h, - h,^^- -i+1- -hj^ 


3(Ui+i-Ui) 3 (u.-u._3_) 

*•' -t* 


h? 

1 


(Ic) 


m 


h . T 
1+1 


i+X-"!" — ‘“i+«i+X> ' 


(Id) 


h . h . u . -u . ^ 

m.= M.+ M. T + 

1 3 1 6 1-1 h. 

1 


(le) 


i 3 ^ ^i ^ ^-+1 — ~ 

1 3 1 6 1+1 h , 

1+1 


(If) 
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2ni. T 4m. u.-u. t 

1-1 . 1^1 1-1 

«i = — - h- 


M.= 

1 


4m, 

1 

1+1 


2m. u . . , -u . 
1+1 , ^ 1+1 1 
- +6 


^i+1 


h? 


i+1 


(ig) 

(Ih) 


Eqs. (lb) or (Ic) lead to a system of N equations for the N+2 

unknowns or m^, respectively. The additional two equations 

are obtained from boundary conditions on m^, %j+l ^+1 

The resulting tridiagonal system for or m^ is diagonally 

3 

dominant and solved by an efficient inversion algorxthm . 

1 2 

Spline 2 for Solving Partial Differential Equations ' 

If the values u^ are not prescribed but represent the solu- 
tion of a quasi-linear second order partial differential equa- 
tion, u =f(u, u , u ), then an approximate solution for u. can 
be obtained by considering the solution of 


This formulation is designated spline 2. If the time derivative 
is discretized in a simple finite-difference fashion, we have 

n+1 n 

u. -u. , 

(l-0)f"+0f"'^-^ ; (2a) 

0=0, explicit; 0=1, implicit; Q=ht Crank-Nicolson . For the 
explicit integration the stability limitations are quite severe, 
see Refs. 1,2 and section VI. Therefore a two-step procedure is 
considered and is given as: 
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1 1 

Step Is = f 


Step 2: 


Example ; 

Consider the linear Burgers equation 


u^+uu^=vu^jj ; u=u(x,t) ; v=v(x,t) 


With (lb) and (Ic) we obtain a system of 3N equations for 
3(N+2) unknowns (see Refs. 1,2 for further details on the 

f 

derivation) . The system (2) can be written as 


A.V^'^^+B.V^‘*‘^+C.V^^^=D.V^+E. [ oV ? T+V^.i ] / 
1 1-1 11 1 1+1 11 1 1-1 1 + 1 -' ' 


where 


A .= 
1 


0 

■1/h. 

3/h? 


0 

0 

1/h. 


Yi 

h^/6 

0 


cc^ 


B 

1 


(1+1/a) 

-3 (l-l/a®)/h® 
1 


a. 


2 ( 1+1/a) /h^ 


Ctg 

(a+l)h^/3 


c.= 

1 


0 

1 


h . , - 
1+1 

-3 


hT 


i+1 


0 

0 

1 


h . - 
1+1 


^+1 

6 

0 


*It is possible to treat the viscous terms (M^) implicitly 
and the convection terms explicitly. As shown in Refs, 1, 
stability of the two-step procedure f'or viscous flows is 

^A number of variations on this system can be derived with 
relations (1) . 


(2b) 


(3a) 


(3b) 


( 0 = 1 ) 

2 , the 
improved . 

the 


9 



and 



L 0 


Pi Ps ■ 

0 0 

0 0 - 

0 6i 

0 0 

0 0 - 


V. = [u., m., M.f , 


(3c) 


^=hi^l/hi ; Yi = Y3 = 6i=0 

a^=l ; ai = 0a^'^^At ; 0 .^= - 0v^'*'^At ; 

pQ=l ; pi= -(l-0)ujAt 7 p2=(l-0)v5At . (3d) 


A significant advantage of the spline 2 formulation is that with 

expressions (1) it is possible to reduce the 3x3 matrix system 

(3) to a scalar set of equations for alone. The details of 

this reduction process are found in Refs. 1,2. 

For equations with two space dimensions such that u^=f(u,u^, 

^y'^xx'^yy^' ^ spline alternating direction implicit (SADI) 

1 2 

procedure has been presented by Rubin and Graves ' . A spline 

successive approximation method can also be simply formulated. 
Both techniques are discussed later in this paper where several 
example problems are presented. 

Truncation Error 

For interior points, the spatial accuracy of the spline 
approximation can be directly estimated from the formulas (lb) 
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and (le) or (If) . Expanding and in Taylor series 

and assuming the necessary continuity of derivatives for u(x,y), 
we obtain, with 

(u ) .=M. + (u^'^) .h?(a®+l)/12(CT+l) 

XX 1 X 'x x' 

(a-1) ( 2 +5 0+2) /ISO 

-(u'^^) ^h^[o®/360+(o-l) ^ (7 ct2_2o+7)/1080] 

+ O(h^) , (4a) 

and 

(Ux) i=nii+(u^'^) ih?o(o-l)/24 + 

t 

+ (u^^)^ht CT[l+a(o-l) ]/180 + 0(h?) . (4b) 

Fyfe^^ has presented similar relations, 

for constant h^, in his collocation analysis of cubic splines 
for the solution of two point boundary value problems. 

Therefore, the spline approximation with a non-uniform mesh 
is second-order accurate for and third-order for m^^. For a 
uniform mesh m^ becomes fourth-order with remaining second- 
order accurate. In the next section a finite-difference expres- 
sion for (u^^) ^ is used to increase the accuracy of and hence 
the overall accuracy of the procedure. With this modification 
this formulation will be termed spline 4. 

+ 

If (Ic) is used to evaluate the truncation error for m^, the 
constant 24 in the second expression on the right-hand side 
becomes 72. For the uniform case, (4b) is recovered in all 
cases. 
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III. SPLINE 4 - DERIVATION AND DISCUSSION 


In order to improve the overall accuracy of the spline 2 
formulation, it is necessary to reduce the order of the trunca- 
tion error for (u ) . in (4a) . Although a number of procedures 

are possible, we have chosen a very simple modification, whereby 

iv . 

the error term in (4a) for (u ) ^ xs approxxmated by a three- 
point discretization for This approximation is first-order 

accurate with a non-uniform mesh and second-order with a uniform 
mesh. Therefore the spline approximation for (u ) . is improved, 
and parallels that for (u ).; i.e., third-order accuracy is 
achieved for a non-uniform grid and fourth-order accuracy for 
uniform mesh. This improvement leads to what is termed spline 
4, or a quintic spline of deficiency three. 

The development of spline 4 is as follows: The expression 

(4a) can be rewritten in the form 

where A= (l+CT®)/a(l+a) ^ . 

The familiar three-point discretization formula is 

^^xx^i= ■STsIt) hi ^^i+1" ] 

-(a-l)h^(M^^^)^/3-h|(l+a=) (M"-^)^/12(l+a)+0(h?) . (5b) 
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Therefore, (4a) or (5a) becomes 


(u^^) ^=M^+ ( A/6 ) ( 1+ a) M^+ 

-7h|(l+cj^) (a-1) (u^) ^/180-hf (u'^^)^Ca®/360 
+ (CT-1) ^ (7a®-2a+7)/1080] + 0(h|) . (5c) 

With (4b) , 

(u^) i=mi+0 (( a-l)h?,ht) ^ (5d) 

and we obtain a uniform higher-order approximation termed 
spline 4. When cp= 1, spline 4 is fourth-order accurate and the 
truncation error of (5b) is smaller than that obtained with 
Hermitian or Fade i;nethods , which are in turn smaller 
than the error obtained with five-point finite-different 
discretizations . 

In the spline 4 procedure the relations (Ib-lh) still apply; 
however, the interpolation polynomial is no longer applicable as 
spline 4 represents a higher-order interpolation. The governing system 


•j* 

It is possible to apply (5b) to (4b) to make (u ) . fourth-order 
even with a non-uniform mesh. 

Higher-order procedures, e.g., spline 6, can be derived in a 
similar manner , and spline 2 is recovered from spline 4 with 
A set equal to zero. 
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remains tridiagonal. Unlike spline 2, where the system can be 
reduced to that for alone, the appearance of off-diagonal 
terms in (5b) restricts the reduction process to a 2x2 
system in (u^,M^) . 

For the linear Burgers equation the system is still of 
the form (3b) with 

Yi= - Vi'^^9AtCTA/6 ; Yg= “ Vi'^^0AtA/6 ; 

ag= - Vi'^^eAt (l-d+a) A/6) ; 

P3= (l-0)vi'^^At(l-(l+a) A/6) ; 

6i= Vi'^^(l-0)At A/6 . (6) 

All other entries in (3c, 3d) are unchanged. 


IV. FINITE-DIFFERENCE THEORY/SPLINE 1 

If the procedures given previously for spline 2 and spline 4 
are repeated for a quadratic polynomial interpolation with both 
derivatives continuous (a quadratic spline of zero deficiency) , 
we find on [x^ I'^i^' 


SA(x)=Ui(x-Xi-i)/h+Ui-i (x^-x)/h+(u^-u^_^-hm^) (x-x^_^) (x^-x)/h®. 


where 


S d ^ • ) ■ » S d X . T ) =u . T , S ' ( X • ) =m . , 
A X X A x-1 x-1 A X 1 


and 


mT = S"Ay~)= - 2(u.-u. T-m.h)/h^ . 
X Ax ' X x-1 1 
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on S^(x.)=u. and S^(Xi)=in. 

we obtain 


4 = s;(xj) = 


- 2 (u.-u. j^T+m.h)/h® 
' 1 1+1 1 


From the continuity of the second-derivative 


and therefore 


wt = m7 

1 1 


■”i= 


The expression for becomes 


M+=M. = (u.^j^-2n.+u._j^)/h= . 


Therefore the quadratic spline of zero deficiency leads to the 
central difference expressions. 


V. STABILITY 


For the linear Burgers equation (3), with u, v constant, the 
interior point stability can be assessed with the von Neumann 
Fourier decomposition of the system (3) for h^=h=constant . 

With V?_j_^=Ti^expI uj(x^+rh) , l=(-l)^, (3) becomes 


T .‘V^ or , 


1 1 


1 1 


1 1 


where G^™T^ is the amplification matrix. The von Neumann 
condition necessary for the suppression of all error growth 


• 15 



requires that the spectral radius p(G^)<l. The eigenvalues of 
G. are , 

X X 

For the one-dimensional equation (3), three numerical pro- 
cedures were considered: (1) convection (m^) and diffusion (M^^) 

explicit, (ii) convection explicit, diffusion implicit (two 
steps required for inviscid stability), and (iii) diffusion and 
convection implicit. With explicit convection, (i) or (ii), 
both divergence and nondivergence forms of the equations have 
been evaluated in Refs. 1 and 2. 

The stability conditions imposed on these schemes is determ- 
ined from 

Uil <1 

(i) Explicit convection and diffusion: 0=0 in (2a, 3). 

Spline | | ®= (l-6@ (1-coscp) (2+coscp) ^)®+c^$^<l, where 

S=vAt/h®, c=uAt/h, $=3sincp/ (2+coscp) , cp=uuh. Necessary stability 
limits are 

(a) 0 < 1/6 , 

(b) c < (3)“'^ , 

_ j, 

(c) R^=c/3=uh/v < 2(3) ^ . (7a) 

These results are more restrictive than the limits found for the 

12 

forward txme central space explxcit f xnite-dxf ference method 
or spline 1, which are 

(a) 8 < 1/2, (b) c < 1 , (c) R^ < 2 . 
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ir 

7 

Spline 4 ; | \ I (l- (5+coscp) (1-coscp) g/ (2+coscp))®+ (Scsincf/ 

(2+coscp))^«l, so that necessary stability limits are 

(a) p<l/4, (b) c<{10)V6 , (c) R^<(40)V3 . (7b) 

Once again these conditions are somewhat more restrictive than 
those obtained with second-order finite-differences (7b) . The 
Fade finite-difference limitation c<(6) VS is even more restric- 
tive, see Appendix of Refs. 1 and 2. It is significant that in 
all cases the explicit method is unconditionally unstable for 
inviscid flow; i.e., p=0. 

(ii) Two-step explicit integration (2b) : 

This procedure, which alleviates the inviscid instability 
found in (i) , is a two-step predictor-corrector method (see 
Refs. 1 and 2) and is similar to the Brailovskaya^^ two-step 
finite-difference technique. For p=0, we obtain 

c < [ (2+coscp) (3sincp) ^ 

This result is more restrictive than the c^l CFL condition found 

for the Brailovskaya finite-difference method. 

For 0 ^ 0 , the effect of diffusion when treated implicitly is 

to improve the inviscid stability limitation. For u-0, the 
. 1 2 

method is unconditionally stable ' . Since the convection terms 

are unchanged, spline 4 has the same stability condition. 

(iii) Implicit convection and diffusion: 

The spline 2 and spline 4 procedures are unconditionally 
stable if 0>l/2 in Eqs. (2), (3). 

17 
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(iv) SADI: 


In Ref. 1, the interior point stability analysis is extended 
to the two-dimensional SADI procedure; unconditional stability 
is demonstrated. 

Although the implicit procedures lead to unconditionally 
stable formulations, as with finite-difference methods, the tri- 
diagonal system may not be diagonally dominant. In this case 

3 

the inversion algorithm may lead to large error growth. Dxag- 
onal dominance can be achieved by a spline adaptation of the 
finite-difference procedure given in Ref. 12. For all the 
problems treated here this modification is unnecessary. In other 
applications it will play a significant role if accurate solu- 
tions are to be obtained. 


VI . RESULTS 

Several model problems have been considered in order to 
evaluate the cubic spline collocation methods presented herein. 

For each of these problems an analytic solution or reliable 
numerical solution is available for comparison purposes. Spline 
interpolation (spline 2 and spline 4) is used to approximate the 
spatial gradients. For the one-dimensional Burgers equation the 
integration procedure outlined in Section II is adopted. Implicit 
or two-step explicit methods are used. For the two-dimensional 
diffusion equation, solutions are obtained with the SADI formula- 
tion. The Laplace equation in cartesian and polar coordinates 
is evaluated with a spline successive approximation procedure. 
Finally, the similarity equations for the flat plate boundary- layer 
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and the two-dimensional stagnation point are solved by direct 
integration of the resulting two-point boundary value problems. 

Solutions are obtained with both uniform and non-uniform 
meshes. Three-point finite-difference calculations are included 
in order to assess the relative increase in accuracy associated 
with the higher-order procedures. The results are presented in 
tabular form so that meaningful comparisons are possible. 

A. Burgers Equation 

The nonlinear Burgers equation (3a), with x=x, u=u(x, t) 
and x=x-(l/2)t, becomes 

u^+(u-l/2)u^=vu^^ (8a) 

with V constant and the boundary conditions 

U”*l as x-» -oo and u-.0 as x-<d . (8b) 

The steady state solution of (8a) is 

u= [l-tanh (x/4v) ]/2 (8c) 

Spline 2 and the finite-difference solutions of (8) have 

* 

been dxscussed in Refs. 1 and 2. Both implicit and two-step 

explicit integration techniques, as outlined in Section II, 

1 2 

have been applied successfully ' . Spline 4 solutions have now 

★ 

been obtained with the implicit and/or two-step procedures of 
(2a), (2b), The system (3) with the coefficients (6) are con- 

sidered. In the actual calculations the 3x3 system (3) is 

* 

The nonlinear coefficient u is treated iteratively or with 
quasi-linearization^» 
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reduced to a 2x2 system. in is eliminated with (le) or 

(If), The "boundary conditions (8b) on u^, specified 

at x=x , with X >3. The boundary conditions on M. are ob- 
max max— i 

tained from the third-order accurate relation 

(u )..,-(u ). =M.,, -M. (9a 

' XX 1+1 XX'^ 1 1+1 1 


where i=l or N. 

The boundary condition (9a) can be applied in two forms. 
These are outlined for the boundary i=l: 

With (u )„ evaluated from (5c), we obtain 

XX 8 

(u^^)i=Mi + (A/6) (M3-(l+a)M3 + CTMi) , (9b) 

where o=h 3 /hg. From the governing equation (8a), 

^^xx^ i/2v-mi/2v 

so that with (If) , 

mi=-h3Mi/3-hgMg/6 + (u^-D/h^ 


and (9b) becomes 

Ml ( 2 v+ CTV A/3+hjg /3 ) +Mg ( h^ /6- v ( 1+ ct) A/3) + ( vA/3 ) M3 -Ug/hg =-l/hg . 

(9c) 

(b) An alternate form of (9a), relating only the two points, 
i=l and i=2, can be derived by evaluating (u )„ from (8a). The 

XX 3 

temporal discretization is given by (2a) . We obtain 

ajM^’^^ + + a3u”'^^ = a^ (9d) 
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where 


3i= At (v+ (Ug+0.5)hg/6 ) ; 
ag=-At(v-(u^-0.25)h3/3) ; 
aa= (l+At (Ug-D/hg) ; 

a4= (Ug+(Ug-l) At/hg) • 

For spline 4 , A= (1+a® ) /a(l+a) ® . For spline 1 , set A=0 so that 
(9c) is second-order accurate. Eq. (9b) is third-order accurate 
for both spline 2 and spline 4. Similar relations are obtained 
for the other boundary, where 

The condition (9c) is independent of the time step At and 
somewhat less cumbersome. It was found that the accuracy of 
the solutions and the time to attain a converged steady state 
solution were virtually insensitive to the choice of the boundary 
condition (9c) or (9d) . This conclusion remains unchanged if 
the higher order effects in (9b), i.e., those terms multiplied 
by A, are treated explicitly in (9c) . In this way (9c) reduces 
to a two-point implicit formula. In several cases the simpler 
spline 2 boundary conditions were applied with the spline 4 
procedure; the solutions always fell between the results of 
spline 2 and spline 4 , but generally closer to those of spline 4. 
Therefore, if simplicity of boundary conditions is desired this 
is a reasonable approximation. 

Typical results are shown, for v=l/8, 1/16, 1/24 on Tables 
1-5. The increase in accuracy as one progresses from the finite- 
difference results to those of spline 2 and finally to spline 4 
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is apparent. This is particularly true with the non-uniform 
meshes of Tables 2 and 3. For the conditions of Table 3, the 
finite-difference calculations with the two-step explicit pro- 
cedure did not converge. An oscillatory behavior was observed 
after 3200 iterations. In certain cases, where h^ is relatively 
large, the nature of the truncation errors (4a, 4b) of spline 2 
and spline 4 is such that a local value obtained with spline 2 
may be as accurate or more accurate than that obtained with 
spline 4. These are exceptional cases, however, and never occur 
for hj^«l. A percentage error plot for the results of Table 1 
is shown on Figure 1. denotes the exact solution (8c). 

Solutions for other v values are of a similar nature and 
therefore have not been included here. 

B. Linear Burgers Equation 

Consider the equation 

u +VU =0 , on 0<x<l , 

with boundary conditions u(l)=l and on x=0, vu +u=0 . The exact 

solution is (x ) =exp ( 1-x ) / v . In some unpublished work by- 

George J. Fix*, it was shown that with this flux boundary 
condition linear finite element theory naturally satisfies the 
required conservation condition at the boundary and therefore 
leads to more accurate solutions than obtained with non-divergence 
versions of spline 2 or conventional finite-difference theory. 

If finite-difference theory is developed in divergence or 
conservation form, the resulting equations are identical with 


*Institute for Computer Applications in Science and Engineering 
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li 

those of the linear, second-order accurate, finite element 
method. If spline 2 is recast in divergence form the solutions 
are considerably more accurate than the non-conservation results 
and also improve upon the conservation finite element (finite- 
difference) calculations. Therefore, with the flux boundary 
condition it appears that divergence form may be required if 
accurate spline solutions are to be obtained. On the other 
hand, if a modified derivative boundary condition was considered 
in lieu of the flux condition, the sensitivity to divergence form 
was no longer apparent. It is possible therefore that the flux 
condition represents a singular case. 

The governing systems of equations and the boundary condi- 
tions for the different formulations are as follows: 

Finite-Dif ference/Non-Diverqence Form 

= 0 (10a) 

At x=l, (lot)) 

At x=0, v(ui_ u_^)/2h + (10c) 

Eliminating u_^ from (10c) with the difference equation (10a), 
we obtain 

2 V® (Ui-u^)/h®+u^ (2 v/h-1) = 0 . (lOd) 

Finite-Difference/Divergence Form = Finite Element 

( vu +u) = 0 

X X 

Therefore, ( i+ig " ^ i-J^ 
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or 


(11a) 






The governing equation (11a) is identical with the non-divergence 
equation (10a) . The alterations appear in the boundary conditions 


At x=l, 

At x=0, (vu +u) i =0. Therefore, 

v(ui~u^)/h + (ui+u^)/2 = 0 . (11b) 

The boundary condition (11b) differs from the non-divergence 
condition (lOd) . 

Spline 2/Non-Diverqence Form 

The governing equation (12a) is combined 

vM^+m^ = 0 , (12a) 

with the spline relations (1) . The boundary conditions are 

Ujj=l, ° * (12b) 

Spline 2/Pivergence Form 

The governing equation (13a) is combined 

(vm+u)^_j_ji= (vm+u)^_j^ (13a) 

with the spline relation (Ic) . The boundary conditions are 



and 

vm^+u^+a( vm^+Uj) = 0 , (13b) 

where a=0 corresponds to the exact boundary value and 

0=1 corresponds to an averaged boundary condition. 
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The results of these calculations are shown on Table 6 




It is seen that the non-conservation (NC) solutions with ten 
mesh points (N=10) are rather poor when compared with either 
the finite element or spline 2 conservation (C) solutions. It 
is significant that the spline divergence solutions, for both 
a=0 and a=l, are considerable improvements over the finite- 
element results. As the number of mesh points increases the 
non-divergence solutions do show some improvement, with spline 2 
more accurate than finite-differences, but these results are 
still less accurate than finite-element solutions. The ten point 
spline 2 divergence form a=l solutions are about as accurate as 
the 50 point finite element results. 

Also shown on the table are ten point solutions with some- 
what modified derivative conditions at x=0. The exact solution 
is unchanged. These derivative boundary conditions were treated 
in much the same manner as the flux condition for each of the 
procedures. For the finite-element solutions an average condition 
was applied. Significantly the large differences between diverg- 
ence solutions no longer occur . The spline solutions are always 
the most accurate, with a small increase in accuracy when diverg- 
ence form is assumed. 

C. Linear Corner Flow 

The two-dimensional diffusion equation 

V ^^'"xx+'"yy^ ' u=u(t,x,y) (14a) 

with the initial condition u(0,x,y)=0 and boundary conditions 
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u(ti>0,0,y^)=l , u(t>0,x>0,0)=l 

u(t,x,y)-+0 as x,y-<o (14b) 

has the exact solution 


u=l-erf X erf Y , (14c) 

where X = ^(R^/t)"^, Y = ^(R^/t)'^ . 

This solution describes the impulsive motion of a right- 
angled corner formed by two infinite flat plates and results of 
the SADI spline 2 calculation have been presented in Refs. 1,2. 

The SADI procedure for the diffusion equation (14a) for both 


spline 2 and spline 4 is given as follows: 

Step 1: = u^ .+ ((u + (u )^.)At/(2R ) 

^ in xn V' XX XT ' yy' 2 .-\/ ' ' e 


(15a) 


Step 2: u^t^ = u^t^ + ((u )^t^ +(u ) At/(2R ) 

^ xj xj ' XX xj ' yY id / “ 


(15b) 


where <«xx> ij= "iJ+< V®> K+i, j) 


and 




and each satisfy a tridiagonal equation of the form 

(lb) . 

^ = (l+o^)/o-y(l+ay) = 


'’x = • hi=Xy-x._y ; k.=y.-y 


The spline 2 formulation is recovered with A =A =0. The boundary 

X y 

conditions for u^ are given by (14b) . The boundary conditions 
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for L. . are obtained from (16) with (u ) . .=(u ) . .=0 on 

13 13 xx' 2.3 yy 

•k 

the boundaries, or from (9a) with (u ^)-,, • obtained from (14a). 

The solution for step 1 is obtained with the tridiagonal 

2x2 system for and u^^ as described by (15a) and (lb). A 

similar procedure for is required for step 2. 

The solution for R =1000 is given on Table 7. A non-uniform 

e 

21x21 mesh with a =c =1,5 was prescribed. The step size At=0.01, 

X y 

The solution is shown for t=2.0. All of the solutions are 
reasonably good for this case, but once again the spline solutions 
are somewhat better. 


D. Laplace Equation 

The Laplace equation 

u +u = 0 ; u=u(x,y) , (17a) 

XX yy 

with the boundary conditions u (0,y) =u (l,y) =0; u (x, 0) =sinTTx; 
lim y -^00 u(x,y)=0 
has the solution 

u(x,y) = (sinnx) exp(-ny) . (17b) 

This boundary layer-like problem was chosen in order to 
evaluate the accuracy of the spline procedures, in particular 
with non-uniform meshes, when large gradients exist only over a 
limited region. In addition* this problem will serve as a 
prototype for spline integration using successive approximation 
(relaxation) procedures. 


This procedure has been demonstrated for the Burgers equation 
by the discussion leading to (9d) . 
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Using the general expression for second derivatives (16 ) , 
Eq. (17a) can be put into a spline form. With the tridiagonal 
relationship for and (lb ) , this leads to a 3x3 system 

for the vector V. . at all interior mesh points: 


A. .V. . t+B. .V. .+C. .V. .^t+D. .V. , .+E..V.JT .= 0 (18) 

x: 1,3-1 13 13 13 1/3+1 13 1-1/3 13 1+1/3 


where 









0 


0 


E 


ij 


0 



X 1 


0 0 
0 1 


The solution is obtained with a successive point relaxation 
procedure. 


V 


(K+1), 

ij 


B 


SC 


A. 

i: i/D-i 


c. + 

i: i/D+i 




^ij^i+1, j- 


(19) 


where the superscript k represents the iteration parameter. 

The system is diagonally dominant and the eigenvalues of the 
amplification matrix, see Section V, are all less than or equal 
to one. The results of this compuation are presented in Table 8, 
The values of u (,5,0) and u(,5,h ) are compared with the exact 

y ® 

solution (17b) , Also included in this table are the results 

obtained with the three-point finite-difference approximation 

for u and u , In order to make a more definitive comparison 
XX yy 

between the spline and finite-difference solutions, the surface 
value of u^ in the latter case was obtained by spline fitting 
the numerical values of u(x,y). In one case noted on Table 8, 
a three-point end difference formula was applied. All of the 
calculations were performed with 10 mesh points in the normal 
or y-direction. In certain cases, spline 2 was used in the y- 
direction and spline 4 in the x-direction. These solutions are 
noted accordingly. 

The spline 4 results are the most accurate in all cases. 

For a uniform mesh the finite-difference and spline 2 results 
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are of equal accuracy as there are no convection effects in the 
problem. Moreover, if the spline 2 and the finite-difference 
solutions are averaged, the spline 4 results are closely approx- 
imated. For a non-uniform mesh the improved accuracy of spline 2 
over the finite-difference approximation is now apparent. 

The spline 4 results are remarkably accurate with o=1.7, 

h =0.1 and v =28.66. For this mesh there are only four points 
s ^ max 

in the region 0<y<l as compared with a uniform mesh (h=0.1) and 
ten points. The coarse mesh, spline 4 results are more accurate 
than the uniform mesh finite-difference solutions. 

The 1.7/. 2 notation for a means that o=1.7 for h^<0.2. 

For h^>0.2, a becomes unity. In this way the mesh width does 
not exceed a specified maximum value. This type of mesh alignment 
is useful in boundary layer problems, where a fine grid is desired 
near the surface, and a uniform but coarser mesh is required in 
the outer inviscid regions. This procedure is also applied for 
the boundary layer solutions in Section VI. F. An error plot is 
given on Figure 2 . 

E. Potential Flow Over a Circular Cylinder 

The governing equation in cylindrical coordinates for the 
potential flow over a circular cylinder is given by 


u + - u + u„^ = 0 

rr r r r^ 99 


( 20 ) 


The boundary conditions are u^(l,9)=0 and lim u(r, 0) -*rcos0. The 

r-»oo 

exact solution of Eq. (20) with these boundary conditions is 
u=(r+^)cos0. Eqs. (17a) and (20) differ only by the appearance 
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of the term which is discretized by the relations (le) or 
(If) . The resulting 3x3 system for and is of the 

form (19) . The coefficient matrices in the present case will 
be somewhat altered by the term. 

The results of the iterative solution are presented in 
Tables 9-11. As in the previous examples, the finite-difference 
solutions are obtained by using three-point central difference 
formulas. In Table 11, the slip velocity on the fore surface 
of the cylinder is presented. The superiority of the spline 
solutions over those resulting from finite-difference discretiza- 
tion is evident. It should be noted that the slip velocity in 
the finite-difference case is obtained by using a three-point 
central difference formula, while the spline solutions require 
only the two-point formula (le) . The higher accuracy of the 
two-point spline formula over the three-point finite-difference 
relations can be of considerable importance for problems with 
derivative boundary conditions. 

F. Similarity Boundary Layers 

The boundary layer equations for the flow 
( p=0) and the two-dimensional stagnation point 
reduced to the following ordinary differential 
appropriate similarity transformations,^'^. 

u "+fu '+p (1-u® ) = 0 

f' = u 

The boundary conditions are 

f(o)=0, u(o)=0, lim u(x)=1.0 

X-*QO 


over a flat plate 
(P=l) can be 
system by using 

(21a) 

(21b) 

(21c) 
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Accurate numerical solutions have been reported in the litera- 

14 

ture [see Rosenhead ], 

In the spline 2 and spline 4 formulation, Eq. (21a) is re- 
duced to a 2x2 system for u^ and and the two-point boundary 
value problem is solved subject to (21c) . For the first-order 
equation (21b), we obtain the following spline approximation 
from (If) : 


h? ^ 

f . , T=f .+h. ^tU.+ (MF. + .5jyiF 

1+1 1 1+1 1 3 1 


i+1' 


( 22 ) 


where MF^=(u')^ for spline 2. For spline 4, the following 
relation to evaluate is easily derived from (le) and (If) : 

MFo=(u')„ • 

Eqs. (22) and (23) give rise to an initial value problem for f^ 
and MF^ which is solved by a marching procedure. Eq. (22) 
leads to third-order accurate expression for f^^^; therefore, for 
non-uniform meshes and third-order accurate solutions, this 
approximation is adequate even for spline 4. For the finite- 
difference solutions, a second-order accurate two-point formula 
for f^, which is consistent with the accuracy of the overall 
scheme, is obtained with the trapezoidal rule. For 3=1, the 
nonlinear term u^ is treated by quasilinearization so that 

k is the iteration parameter. 
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The results of these computations for both uniform as 
well as non-uniform meshes are tabulated in Tables 12-16, The 
shear at the wall is proportional to f"(o) and this term has 
been evaluated for a variety of meshes. The results are given 
on Tables 15 and 16 for p=0 and 0=1, respectively. Finite- 
difference solutions for u^ are obtained by using the three- 
point central difference approximation. 

As noted previously, the notation o=1.8/l means that o=1.8 
until h. reaches 1.0, at which point o=1.0. h_ is the first 

X “ 

mesh width off the wall x=0; N is the total number of mesh 
points. Ng is the number of mesh points in the boundary layer 
defined by x<6 . At x=6, |u-1.0|<10 . 

B=0 Blasius solution ; 

The spline 4 solution for N=61, hg=0.1 and o=1.0 is almost 
identical with the "exact" solution of f " (o)=0. 469600. If 
spline 2 boundary conditions are used with a spline 4 interior 
point formulation,f"(o)=0.46 9608. As previously noted, this value 
lies between the spline 2 and spline 4 results. With cp= 1.8/2, 
ha=0.5, N=21 and only 5 points in the boundary layer (Ng^S) , 
the spline 2 value of f"(o) is in error by only 2%. For the 
larger hg values the spline 2 solutions are even more accurate 
than those found with spline 4. Similar behavior was observed 
with Burgers equation in Section VI .A, An error plot is given 
on Figure 3. 

B=1 stagnation point flow ; 

For 3=0, the exact solution has u "'(o) =u^^ (o) =0 and there- 
fore the inherent lower-order accuracy of the finite-difference 
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calculation is somewhat obscured near the wall x=0. For the 
stagnation point solution where f"(o) =1,232588, the improvement 
associated with the spline formulation is clearly demonstrated. 
Therefore, it would appear that spline integration should be 
extremely useful for boundary layer problems. 


G. Non-Similar Boiuidary Layer Analysis 

As a final test of the spline procedures the constant 
pressure boundary layer equations written in physical variables 
(x,y) were considered: 


uu + vu 
X y 



u 


yy 


u + 

X 


V 


y 


0 


The boundary conditions are 


y=0 : u = V = 0 
y»l u -* 1 


The initial conditions were given by 

u(o,y) = 1.0, y^Q and u(o,o)=0 
v(o,y) = 0 

The equations were integrated for a Reynolds number Re=10® and 
a non-uniform mesh of ten points normal to the surface. The 
solution for the normalized skin friction is shown on Figure 4, 
The value Ns denotes the actual number of points within the 
boundary layer. The same criteria of Section VI. F was applied. 
As the boundary layer grows with distance x, Ng increases. 

With 6 to 7 points in the final boundary layer profiles, the 
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spline 4 solutions are quite accurate 


VII . SUMMARY 

It has been demonstrated that higher-order calculation 
procedures using cubic spline collocation provide accurate 
solutions to a number of model problems. The spline methods 
termed spline 2 and spline 4 can be used for two-point boundary 
value problems, as well as implicit, explicit, two-step, ADI 
and iterative integration procedures. 

Spline 4 is fourth-order accurate with a uniform mesh and 
third-order with a moderate non-uniform mesh. Spline 2 is 
second-order accurate for diffusion terms and fourth-order 
(third-order) for convection with a uniform (non-uniform) mesh. 
Derivative boundary values are obtained directly without the 
need for end differencing. For implicit linear systems, the 
spline methods remain unconditionally stable. 

The results confirm the higher-order accuracy of the spline 
methods and lead to the hopeful conclusion that accurate solu- 
tions for more practical flow problems can be obtained with 
relatively coarse non-uniform meshes. 


There has been no attempt to optimize the temporal inte- 
gration procedure so as to minimize computer times or increase 
temporal accuracy. The finite-difference calculations run 20% 
to 25% faster than the spline integrations. When spline fitting 
is used to evaluate finite-difference derivatives, as in 
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Section VI. C, the computer times are comparable. It is antici- 
pated that the reduced mesh requirements with these spline 
methods will result in a net improvement in computer storage 
and time. 
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TABLE I 


SOLUTION OF BURGERS EQUATION 
z/ = l/8, O- = l.0, 31 EQUALLY 
SPACED POINTS 


X 

FO. 

SPLINE 2 

SPLINE 4 

EXACT 

0 

0.5000 

0.5000 

0.5000 

0.5000 

-0.2 

0.6999 

0.6860 

0.6900 

06900 

-0.4 

0.8447 

0.8290 

0.8322 

0.8320 

-0.6 

0.9269 

0.9160 

0.9170 

0.9170 

-0.8 

0.9673 

0.9620 

0.9609 

0.9610 

- 1.0 

0.9857 

0.9830 

0.9820 

09820 

- 1.2 

0.9938 

0.9930 

0.9918 

0.9920 

- 1.4 

0.9973 

0.9970 

0.9963 

0.9960 

- 1.6 

0.9988 

0.9990 

0.9983 

0.9980 

- 1.8 

0.9995 

0.9990 

0.9993 

0.9990 
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TABLE 2 

SOLUTION OF BURGERS EQUATION 
l/ = 1/8,0- = 1.2, 15 POINTS 


X 

FO. 

SPLINE 2 

SPLINE 4 

EXACT 

0 

0.5000 

0.5000 

0.5000 

0.5000 

-0.3859 

0.9510 

0.8214 

0.8297 

0.8240 

^0.8494 

1 .0030 

0.9778 

0.9654 

0.9676 

-1 .4060 

0.9990 

1 .0004 

0.995 1 

0.9964 

-2.0750 

1.0 

1.0 

0.9989 

0.9997 

-2.8770 

1.0 

1.0 

0.9999 

1.0 

-3.8420 

1.0 

1.0 

0.9996 

1.0 

-5.000 

1.0 

1.0 

1.0 

1.0 
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TABLE 3 


SOLUTION OF BURGERS 

equation: 1/ =1/8,0- = 1.8, 

15 POINTS 


X 

SPLINE 2 

SPLINE 4 

EXACT 

0 

0.5000 

0.5000 

0.5000 

-0.0662 

0.5740 

0.5691 

0.5659 

-0.1855 

0.6986 

0.6858 

0.6774 

-0.4001 

0.8695 

0.8452 

0.8321 

-0.7864 

1.0012 

0.9689 

0.9587 

-1 .4818 

1.0165 

1 .0083 

0.9973 

-2.7334 

1 .0267 

1.0257 

1.0 

-4.9864 

1.0 

1.0 

1.0 
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TABLE 4 SOLUTION OF BURGERS 

equation: v=i/i 6, <r=i.o, 

19 EQUALLY SFV\CED POINTS 



FO. 

SPUNE 2 

SPLINE 4 

EXACT 

0 

0.5000 

0.5000 

0.5000 

0.5000 

-0.2 

0.9000 

0.8231 

0.8356 

0.8320 

-0.4 

0.9878 

0.9641 

0.9617 

0.9608 

-0.6 

0.9986 

0.9952 

0.9916 

0.9918 

-0.8 

0.9998 

0.9995 

0.9982 

0.9983 

- 1.0 

1.0 

0.9999 

0.9996 

0.9997 

- 1.2 

1.0 

1.0 

0.9999 

0.9999 
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TABLE 5 SOLUTION OF BURGERS 

equation: V = 1/24, <r= 1.2, 
31 POINTS 



SPLINE 2 

SPLINE 4 

EXACT 

0 

0.5000 

0.5000 

0.5000 

-0.0688 

0.6936 

0.6957 

0.6955 

-0.1514 

0.8618 

0.8606 

0.8602 

-0.2505 

0.9586 

0.9526 

0.9529 

-0.3695 

0.9928 

0.9876 

0.9883 

-0.5122 

0.9995 

0.9975 

0.9979 

-0.6835 

1.0 

0.9996 

0.9997 

-0.8890 

1.0 

1.0 

1.0 

-1.1356 

1.0 

1.0 

1.0 

-4.9582 

1.0 

1.0 

1.0 
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TABLES LINEAR BURGERS EQUATION 


N 

FINITE-ELEMENT 

SPLINE 

FINITE -DIFFERENCE 

iO 

165.38 

37378(NC) 

14.7 

25 

150.923 

100.632 (NO 


50 

149.034 

132.668 (NO 


100 

148.568 

144.137 (NO 

136.022 

10 


138. 6 (C) 

— 

10 


149.1 15 (C AVG.) 

— 


(I'Ux +2U) 

• e''" AT THE BOUNDARY 

10 

165.38 

151.552 (NO 
148.486 (C) 

158.969 


(vUx) • - e 

AT THE 80UN0ARY 

10 

I65i38 

145.442 (NO 
148.342 (0 

139.296 


EXACT SOLUTION Uq • 148.4 , v = 1/5 



TABLE 7 LINEARIZED CORNER FLOW 


METHOD 

Y 

0.0025 

0.0081 

0.0493 

0.2575 

I.3II7 

6.6485 

EXACT SOLUTION 

<r »l.5 , At » 0.01 
k2 = h2=-00l 

Nx ■ Ny *21 
t =2.0 

0.0025 

0.0081 

0.0493 

0.2575 

I.3II7 

6.6485 

0.9990 

0.9968 

0.9822 

0.9685 

0.9685 

0.9685 

0.9968 

0.9896 

0.9424 

0.8978 

0.8978 

0.8978 

0.9822 
0.9424 
0.6820 
0.4361 
0.4361 
0.436 1 

0.9685 

0.8978 

0.4361 

0.0001 

0.0 

0.0 

0.9685 

0.8978 

0.4361 

0.0 

0.0 

0.0 

0.9685 

0.8978 

0.4361 

0.0 

0.0 

0.0 

FINITE DIFFERENCE 

0.0025 

0.0081 

0.0493 

0.2575 

I.3II7 

6.6485 

0.9990 

0.9967 

0.9818 

0.9681 

0.9681 

0.9681 

0.9967 

0.9893 

0.9408 

0.8966 

0.8966 

0.8966 

0.9818 

0.9408 

0.6724 

0.4275 

0.4273 

0.4273 

0.9682 

0.8967 

0.428 1 

0.0007 

0.0004 

10.0004 

0.9681 

0.8967 

0.4279 

0.0004 

0.0 

0.0 

0.9681 

0.8967 

0.4279 

0.0004 

0.0 

0.0 

SPLINE 2 

0.0025 

0.0081 

0.0493 

0.2575 

I.3II7 

6.6485 

0.9990 

0.9968 

0.9821 

0.9683 

0.9683 

0.9683 

0.9968 

0.9895 

0.9418 

0.8973 

0.8973 

0.8973 

0.9821 

0.9419 

0.6790 

0.4330 

0.4332 

0.4332 

0.9684 
0.8974 
0.4336 
-.0005 
i - .0003 
-.0003 

0.9684 

0.8974 

0.4338 

-.0003 

0.0 

0.0 

0.9684 

0.8974 

0.4338 

-.0003 

0.0 

0.0 

SPLINE 4 

0.0025 

0.0081 

0.0493 

0.2575 

I.3II7 

6.6485 

0.9990 

0.9968 

0.9822 

0.9865 

0.9684 

0.9684 

0.9968 

0.9895 

0.9423 

10.8977 

10.8977 
iO.8977 

0.9822 

0.9424 

0.6820 

0.4360 

0.4358 

0.4358 

0.9685 

0.8979 

10.4366 

j0.0007 

0.0003 

,10.0003 

ii 

0.9685 

0.8978 

0.4364 

0.0003 

0.0 

0.0 

0.9685 

0.8978 

0.4364 

0.0003 

0.0 

0.0 




TABLE 8 

SOLUTION OF THE LAPLACE 
EQUATION 


METHOD 

-UJIS.0) 

U (.Ahj) 

O’ 

kg 

hg 

^MAX 

EXACT SOLUTION 

3.142 

0.7304 

1.0 

0.1 

0.1 

1.0 

FINITE DIFFERENCE 

3.123 

0.7322 





SPLINE 2 

3.164 

0.7286 





SPLINE 4X, SPLINE 2Y 

3.154 

0.7295 





SPLINE 4 

3.142 

0.7304 





EXACT SOLUTION 

3.142 

0.5335 

1.0 

0.1 

0.2 

2.0 

SPLINE 4 

3.162 

0.5329 





SPLINE 4X, SPUNE 2Y 

3. 193 

0.5277 





FINITE DIFFERENCE 

2.828* 

0.5401 





EXACT SOLUTION 

3.142 

0.9691 

1.7 

0.1 

0.01 

2.86 

FINITE DIFFERENCE 

3. 189 

0.9686 





SPLINE 2 

3. 175 

0.9687 





SPLINE 4X, SPLINE 2Y 

3. 162 

0.9689 





SPLINE 4 

3. 137 

0.9691 





EXACT SOLUTION 

3. 142 

0.9391 

1.7 

0.1 

0.02 

2.73 

FINITE DIFFERENCE 

3. 193 

0.938 1 





SPLINE 2 

3. 177 

0.9384 





SPLINE 4X, SPLINE 2Y 

3. 164 

0.9387 





EXACT SOLUTION 

3. 142 

0.7304 

1.7 

0.1 

0.1 

28.66 

FINITE DIFFERENCE 

3.220 

0.7233 





SPLINE 4X, SPLINE 2Y 

3.185 

0.7263 





SPLINE 4 

3.130 

0.7313 





EXACT SOLUTION 

3.1416 

0.9691 

17/2 

0.1 

0.01 

1.063 

FINITE DIFFERENCE 

3.1694 

0.9688 





SPLINE 2 

3.1676 

0.9688 





SPLINE 4X, SPLINE 2Y 

3.1551 

0.9689 





SPLINE 4 

3. 1404 

1 

0.9691 


1 




* EVALUATED BY 3 -POINT END- DIFFERENCE FORMULA 
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TABLE 9 POTENTIAL FLOW OVER A CIRCULAR 

CYLINDER 


METHOD 


0 

TT/IO 

2tt/I0 

3tt/I0 

4 tt/IO 

EXACT SOLUTION 

1.0500 

2.0024 

1.9044 

1.6200 

1.1770 

0.6188 

^9 = tt/IO 

1.2795 

2.0611 

1.9602 

1.6674 

I.2II5 

0.6369 

ha =0.05 

1.9428 

2.4575 

2.3372 

1.9881 

1.4445 

0.7594 

O- =1.7 

3.8596 

4.1187 

3.9171 

3.3321 

2.4209 

1 .2727 

'5-3285 

9.3991 

9.5055 

9,0403 

7.6901 

5.5872 

2.9374 

FINITE DIFFERENCE 

1.0500 

1 .8626 

1.7714 

1 .5069 

10948 

0.5757 


1.2795 

1 .9365 

1.8417 

1.5667 

1.1383 

0.5985 


1.9428 

2.3597 

2.2442 

1.9091 

1.3871 

0.7293 


3. 8 596 

4.0573 

3.8587 

3.2824 

2.3848 

1.2538 


9.3991 

9.4618 

8.9987 

7.6547 

5.5615 

2.9239 

SPLINE 2 

1.0500 

1.9249 

1.8307 

1.5573 

1.1314 

0.5948 ! 


1.2795 

1.9821 

1.8851 

1.6035 

1.1650 

0.6125 

' 

1.9428 

2.3726 

2.2565 

1.9195 

1.3946 

0.7332 

; 

3.8596 

4 .0285 

3.8314 

3.2591 

2.3679 

1 .2449 


9.3991 

9.4254 

8.9641 

7.6253 

5.5401 

2.9126 

SPLINE 4 

1.0500 

2.0089 

1.9106 

1.6253 

1.1808 

0.6208 

1 

1.2795 

2.0677 

1.9665 

1.6728 

1.2154 

0.6390 

1 

1.9428 

2 .4 639 

2.3433 

1.9934 

1.4483 

0.7614 


3.8596 

4.1185 

3.9169 

3.3319 

2.4208 

1 .2727 

:i 

9.3991 

9 .4760 

9.0122 

7.6663 

5.5699 

2 .9283 


TABLE 10 POTENTIAL FLOW OVER A CIRCULAR 

CYLINDER 


METHOD 


0 

ir/10 

27T/IO 

3ir/IO 

I 

4ir/IO 

EXACT SOLUTION 

1.0500 

2.0024 

1.9044 

1.6200 

1.1770 

0.6188 

l9 = tt/20 

1.279 5 

2.0611 

1.9602 

1.6674 

1.2115 

0.6369 

ho *0.05 

1.9428 

2.4575 

2.3372 

1.9881 

1.4445 

0.7594 

! <T *1.7 

3.8596 

4.1187 

3.9171 

3.3321 

2.4209 

1 .2727 

*^MAX “ 

9.3991 

9.5055 9.0403 

7.6901 

5.5872 

2.9374 

FINITE DIFFERENCE 

1.0500 

1 .9246 

1.8304 

1.5571 

1.1313 

0.5949 


1.2795 

1.9830 

1.8859 

1.6043 

1.1657 

0.6130 ! 


1.9428 

2.3816 

2.2650 

1.9268 

1 .4000 

0.7361 1 


3.8596 

4.0551 

3.8566 

3.2807 

2.3836 

1 .2532 


9.3991 

9.4506 

8.9881 

7.6457 

5.5549 

2.9204 

SPLINE 2 

1.0500 

1.9413 

1.8463 

1.5706 

I.I4I2 

0.6001 


1.2795 

1.9986 

1.9008 

1.6169 

1.1749 

0.6178 


1.9428 

2.3898 

2.2728 

1.9334 

1.4048 

0.7387 


3.8596 

4.0474 

3.8494 

3.2745 

2.3791 

1 .2509 


9.3991 

9.4399 

8.9778 

7.6370 

5.5486 

2.9171 

SPLINE 4 

1.0500 

2.0094 

I.9III 

1.6257 

1.1813 

0.6212 


1.27 95 

2.0682 

1 .9670 

1.6732 

1.2158 

0.6394 


1.9428 

2.4643 

2.3437 

1.9938 

1.4487 

0.7618 


3.8596 

4.1188 

3.9172 

3.3323 

2.4211 

1.2730 


9.399 1 

9.4762 

9.0 124 

7.6664 

5.5700 

2.9284 






TABLE II 


SLIP VELOCITY ON 
THE FRONT OF A 
CIRCULAR CYLINDER 

7t/IO 


9 

F.O. 

SPLINE 2 

SPLINE 4 

EXACT 

tt/IO 

-0.564 

-0.594 

-0.620 

-0.628 

2tt/IO 

-1.073 

- 1.130 

-1 .179 

-1.176 

3ir/IO 

-1.477 

-1.555 

- 1 .623 

-1.618 

4ir/IO 

-1.736 

-1.628 

- 1 .908 

-1.902 

5tt/IO 

-1.825 

- 1.922 

-2.006 

-2.000 
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TABLE 12 


BLASIUS profile: <r=IO, h2=0.l, N=6I 



f 


X 

F.D. 

SPLINE 2 

SPLINE 4 

EXACT 

FO. 

SPLINE 2 

SPLINE 4 

EXACT 

0.1 

0.002348 

0.002348 

0.002348b.002348 

0.046967 

0.046962 

0.046959 

0.046959 

0.2 

0.009392 

0.009392 

0.009391 

0.009391 

0.093923 

0.093908 

0.093905 

0.093905 

0.4 

0.037555 

0.037551 

0.037549 

0.037549 

0. 187648 

0.187604 

0. 187605 

0.187605 

0.6 

0.084399 

0.084387 

0.084386 

0.084386 

0.280651 

0.280563 

0.280576 

0.280575 

0.8 

0.149697 

0.149673 

0.149675 

0.149675 

0.372076 

0.371934 

0.371964 

0.371963 

1.0 

0.233026 

0.232982 

0.232990 

0.232990 

0.460788 

0.460583 

0.460633 

0.460632 

1.5 

0.5151 II 

0.514990 

0.515032 

0.51503 1 

0.661735 

0.661379 

0.661474 

0.661473 

o 

cvi 

0.886938 

0.886707 

0.886798 

0.886796 

0.817023 

0.816600 

0.816695 

0.816694 

4.0 

2.784256 

2.783770 

2.783890 

2.783885 

0.997824 

0.997790 

0.997771 

0.997770 

6.0 

4.783607 

4.783110 

4.783220 

4.783217 

1.0 

1.0 

1.0 

1.0 



TABLE 13 


BLASIUS PROFILE ! a =1.0, hg = 1.0, N = 21 



f 

f' 

X 

FD. 

SPLINE 2 

SPLINE 4 

EXACT 

F.O. 

SPLINE 2 

SPLINE 4 

EXACT 

1.0 

0.23859 

0.23768 

0.23490 

0.23299 

0.47718 

0.45853 

0.46175 

0.46063 

2.0 

0.9035 1 

0.8918 1 

0.8883 1 

0.88679 

0.85265 

0.81125 

0.81795 

0.81669 

3.0 

1 .82705 

1.79403 

1 .80274 

1.79557 

0.99444 

097059 

0.96996 

0.96905 

4.0 

2.82470 

2.77982 

2.78658 

2.78389 

1 .00085 

0.99948 

09970 1 

0.99777 

5.0 

3.82500 

3.77942 

3.78874 

3.78323 

0.99975 

099999 

1.00001 

0.99994 

6.0 

4.82492 

4.77947 

4.78555 

4.78322 

1 .00010 

1.00002 

0.99989 

1.0 

20.0 

18.8249 

■ 

18.7795 

18.7857 

18.78322 

1.0 

1.0 

l.O 

1.0 



TABLE 14 


BLASIUS profile: cr=l.8/| ,h2=0.5, N=2I 



f 

f' 

X 

FO. 

SPLINE 2 

SPLINE 4 

EXACT 

F.O. 

SPLINE 2 

SPLINE 4: 

EXACT 

0.5 

0.05910 

0.05835 

0.05888 

0.05864 

0.23642 

0.23230 

0.23477 

0.23423 

1.4 

0.45557 

0.45051 

0.45722 

0.45072 

0.64461 

0.61651 

0.62568 

0.62439 

2.4 

1 .2455 1 

1 .2253 1 

1 .24067 

1.23153 

0.93527 

0.89654 

0.90300 

0.90107 

3.4 

2.21456 

2.17478 

2.19064 

2.18747 

1 .00284 

0.99037 

0.98802 

0.98797 

4.4 

3.21568 

3.16972 

3.18422 

3.18338 

0.99940 

1 .00010 

0.99898 

0.99940 

5.4 

4.21548 

4. 16985 

4. 18386 

4.18322 

1.00020 

1. 00001 

1.00003 

0.99999 

19.4 

18.2155 

18.1698 

18. 18380 

18.18322 

1.0 

1.0 

1.0 

1.0 



TABLE 15 


f"(0) FOR BLASIUS EQUATION 


Xmax 

\)2 

0* 

FINITE 

DIFFERENCE 

SPLINE 2 

SPUNE 4 

Ng/N 

6.0 

0.1 

1.0 

0.4697265 

0.469634 

0.469601 

61/61 

20.0 

1.0 

1.0 

0.528041 

0.475357 

0.476359 

7/21 

5.6665 

0.05 

1.5 

0.516646 

0.470718 

0.466048 

11/ II 

11.3330 

6.4344 

0.1 

0.2 

1.5 

1.5 

0.6049558 

0.498214 


0.493598 

0.455623 

9/11 

7/8 

13. 365 
16.063 

0.01 

0.05 

1.8/, 

1.8/, 

0.474643 

0.473974 


0.469188 

0.469438 

13/21 

10/21 

19.400 

0.5 

1.8/1 

0.479715 

0.466839 

0.469509 

7/21 

37.020 

0.5 

1.8/2 

0.551803 

0.460823 

0.477930 

5/21 

53.936 

0.5 

1.8/3 

0.827648 

0.506798 

0.523256 

5/21 


f"(0) *0.469600 (ROSENHEAo'*®’) 
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TABLE 16 


STAGNATION POINT FLOW 

(a)f"{0) 


Xmax 

hz 

<7 

FINITE 

DIFFERENCE 

SPLINE 2 

SPLINE 4 

N 5 /N 

6.0 

0.1 

1.0 

1.23257 

1.23227 

1.23258 

51/61 

20.0 

1.0 

1.0 

1.07167 

1.20612 

1.20882 

6/21 

9.4448 

0.001 

1.8/, 

1.26353 

1.23604 

1.23299 

16/21 

19.40 

0.5 

1.8/, 

1.24031 

1.22764 

1.23617 

6/21 


f"(0) •1.232588 ( ROSENHEAD''®') 


(b) f (hz) 


Y 

'^MAX 

hz 

<T 

FINITE 

DIFFERENCE 

SPLINE 2 

SPLINE 4 

' 

ROSENHEAD*'®' 

6.0 

0.1 

1.0 

0.005915 

0.005995 

0.005996 

0.005996 

20.0 

1.0 

1.0 

0.390440 

0.436393 

0.450482 

0.459227 

19.40 

0.5 



'•®/l. 

0.128780 

0.132622 

0.135410 

0.133585 




o 



FIG. I NON-LINEAR BURGERS 


BURGERS EQUATION 

O FINITE DIFFERENCE 
({) SPLINE 2 
• SPLINE 4 


O 


O 


O 

i> ^ 



0.8 1.0 1.2 1.4 1.6 


equation: v«I/8, ct-I 








o (76.24) 


0(28.82) 


0(17.5) 
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BLASIUS EQUATION 

O FINITE DIFFERENCE 
4) SPLINE 2 
• SPLINE 4 


12 


O 


10 


o 

= 01 


8 




o 

o 


2h 






4> 


6 

—f— 


rO(. 02693) 


5/21 5/21 7/8 7/21 

1.8/2 1.8/3 1.5 1.0 


"T I T 

7/21 9/11 10/21 ll/n 13/21 61/61 N 

1.8/1 1.5 1.8/1 1.5 1.8/1 1 


6 


FIG. 3 ERROR PLOT ! BLASIUS EQUATION 
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NASA- Langley, 1976 CR— 2653 



FINITE DIFFERENCE 

SPLINE 4 

O’ =1.5 
hg = .001 

N =10 
Rg = IO® 


0.8 1.0 1.2 1.4 1.6 X 

7 Ne 

BOUNDARY LAYER SOLUTION - PHYSICAL 



